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AN ORDINARY DIFFERENTIAL EQUATION MODEL 
FOR FISH SCHOOLING 


Takeshi UcHiTANEt, Ta Viet Ton^, Atsushi Yagi^* 
Osaka University 

This paper presents a stochastic differential equation model for 
describing the process of fish schooling. The model equation always 
possesses a unique local solution, but global existence can be shown 
only in some particular cases. Some numerical examples show that 
the global existence may fail in general. 


1. Introduction. We are interested in describing the process of fish 
schooling by the ordinary differential equations. A model written in terms 
of ODE is very useful. First, the rules of behavior of individual animals can 
be described precisely. Second, many techniques which have been developed 
in the theory of ODE can directly be available to analyse their solutions 
including asymptotic behavior and numerical computations. 

We will regard the fish as particles in the space The direction in which 
a fish proceeds is regarded as its forward direction. As for the assumptions 
of modeling, we will follow the idea presented by Camazine-Deneubourg- 
Franks-Sneyd-Theraulaz-Bonabeau [3] which is also based on empirical re¬ 
sults Aoki [1], Huth-Wissel [6] and Warburton-Lazarus [11]. In the mono¬ 
graph [3, Chapter 11], they have made the following assumptions: 

1. The school has no leaders and each fish follows the same behavioral 
rules. 

2. To decide where to move, each fish uses some form of weighted average 
of the position and orientation of its nearest neighbors. 

3. There is a degree of uncertainty in the individual’s behavior that re¬ 
flects both the imperfect information-gathering ability of a hsh and 
the imperfect execution of the fish’s actions. 

We remark that similar assumptions, but deterministic ones, were also in¬ 
troduced by Reynolds [9]. 
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As seen in Section 2, we formulate the motion of each individual by a 
system of deterministic and stochastic differential equations. The weight of 
average is taken analogously to the law of gravitation. That is, for the i-th 
fish at position Xj, the interacting force with the j-th one at Xj {i ^ j) is 
given by 


—a 
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{\\xi-Xj\\/r)P 


\\Xi 


1 

Xj||/r)T 


{Xi 


Xj), 


where 1 < p < q < oo are some fixed exponents and r > 0 is a critical radius. 
This means that if Xi and Xj are far enough that \\xi — Xj\\ > r, then the 
interaction is attractive; conversely, if it is opposite \\xi — Xj\\ < r, then the 
interaction is repulsive. The exponents p, q and the radius r may depend on 
the species of animal. The larger p and q are, the shorter the relative range 
of interactions between two individuals. 

A similar weight of average is used for the orientation matching, too, i.e.. 
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Here, Vi and Vj denote velocities of the z-th and j-ih. animals, respectively. 

Several kinds of mathematical models have already been presented, in¬ 
cluding difference or differential models. Vicsek et al. [10] introduced a sim¬ 
ple difference model, assuming that each particle is driven with a constant 
absolute velocity and the average direction of motion of the particles in its 
neighborhood together with some random perturbation. Oboshi et al. [7] 
presented another difference model in which an individual selects one ba¬ 
sic behavioral pattern from four based on the distance between it and its 
nearest neighbor. Finally, Olfati-Saber [8] and D'Orsogna et al. [4] con¬ 
structed a deterministic differential model using a generalized Morse and 
attractive/repulsive potential functions, respectively. 

In this paper, after introducing the model equations, we shall prove local 
existence of solutions and in some particular cases global existence, too. We 
shall also present some numerical examples which show robustness of the be¬ 
havioral rules introduced in [3, Chapter 11] for forming a swarm against the 
uncertainty of individual’s information processing and executing its actions. 

In the forthcoming paper, we are going to construct a particle swarm 
optimization scheme on the basis of the behavioral rules of swarming animals 
which can spontaneously and successfully find their feeding stations. 

The organization of the present paper is as follows. In the next section, 
we show our model equations. Section 3 is devoted to proving local exis¬ 
tence of solutions. Section 4 gives global existence for both deterministic 
and stochastic cases but the number of animal is only two. Some numerical 
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examples that suggest global existence is not true in general are presented 
in Section 5. 

2. Model Equations. We consider motion of N fish. They are re¬ 
garded as moving particles in the space {d = 1,2,3,...). The position 
of the 1-th particle is denoted by x* = Xi{t) (i = 1,2,..., N). Its velocity is 
denoted by Vi = Vi{t) (1 = 1,2,..., N). Our model is then given by 

dxi = Vidt + (Tidwi, 



( 2 . 1 ) 



V. 


The first equation is a stochastic equation on Xi, where aidwi{t) denotes a 
noise resulting from the imperfectness of information-gathering and action of 
the fish. In fact, {wi{i),t > 0}(1 = 1,..., N) are independent d - dimensional 
Brownian motions defined on a complete probability space with filtration 
{ld,F, {J-t}t>o,P) satisfying the usual conditions. The second one is a deter¬ 
ministic equation on Vi, where 1 < p < q < <x> are hxed exponents, r > 0 is a 
fixed radius and a, fd are positive constants. Finally, Fi{t,Xi,Vi) denotes an 
external force at time t which is a given function defined for [xi, Vi) with val¬ 
ues in M'^. It is assumed that Fi{t,Xi,Vi) {i = 1,... ,N) are locally Lipschitz 
continuous. 

In what follows, for simplicity, we shall put ai = ar^, jdi = /3r^, 7 = 

Then, the system (2.1) is rewritten in the form 


( 2 . 2 ) 



for 1 = 1 ,..., W 
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3. Local Solution. We set the phase space 

R{N) = {(xi,... ,XN,vi, ... , vn) G X I Xi 7^ Xj 
(1 <ij < N,i^j)}. 

Since all the functions in the right hand side of (2.2) are locally Lipschitz 
continuous in the existence and uniqueness of local solutions to (2.2) 

starting from points belonging to this phase space are obvious in both de¬ 
terministic and stochastic cases, see for instance [2, 5]. Thus, we have 

Theorem 3.1. For any initial value 

(xi(0),... ,XAr(0),ui(0),... ,UAr(0)) G M(iV), 

(2.2) has a unique local solution defined on an interval [0, r) with values in 
R{N), where r < 00 and if t < 00 it is an explosion time. 

4. Global solution in some particular cases. In this section, we 
shall consider the case where N = 2 and prove global existence for (2.2). 
First, the deterministic case (i.e., ui = (T 2 = 0) is treated with null external 
forces Fi = F 2 = 0. Second, the stochastic case (i.e., ai + a 2 > 0) is treated 
but under the restriction that d and q satisfy the relations d > max {(7 —4, 2} 
and q > 2 (therefore, in particular, d > 2 ). 

4.1. Deterministic case (ui = cr 2 = 0) 

(4.1) 

dxi 

lit 

dvj _ ai{xi-Xj) ^ aij{xi-Xj) 
dt llxj — Xjlp ||xj — Xjll'i 

where i,j = l, 2 ,i ^ j. 

Theorem 4.1. Let 1 < p < q < 00 and q > 2 . Then, for any initial value 
(x°,u°) G M(2), (4.1) has a unique global solution {x{t),v{t)) with values in 
M(2). 

Proof. As stated in Theorem 3.1, there is a unique solution { xi { t ), X2 { t ), 
vi{t),V2{t)) to (4.1) defined on an interval [0,ri), where ri denotes the ex- 


The system (2.2) has the form 


fiijvi - Vj) _ iiljvi - Vj) 

I I rp . rp . I IP II rp . rp . | | O' 

I |iL2 J I I I \'^'i '^3 I I 
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plosion time. On [0, ti), (4.1) is equivalent to 
d{xi + X 2 ) 

dt 

d{vi + V 2 ) _ 
dt ~ ’ 
d(xi - X2) 

d{vi - V2) _ 
dt 

-2 

Thus, 


Oil 


ai7 


\\xi-X 2 \\P IIX1-X2II'? 

/ 3 i , / 3 i 7 


\\XI-X2\\P \\xi-X2\\^ 


(xi - X2) 
(ui - V2). 


(4.2) 


' Xl{t) + X 2 {t) 
Vl{t) + V 2 {t) 
d{xi - X2) 

dt 

d{vi - V2) 
dt 


[ui( 0 ) +?; 2 ( 0 )]t+ a;i( 0 ) + 3:2(0), 

t’l(O) + U 2 ( 0 ), 

Vl - V2, 

_2 «! _ Q:i7 

[llxi - X2\\P \\xi - X2||^ 

- 2 r _L ^7 

l_||xi-X 2 ||P ||xi-X 2 ||‘' 


(xi - X2) 
(ui - V2). 


So we put C = xi — X2 and rj = vi — V2- In order to prove that ri = 00 , it 
suffices to show that the solution starting in = {^ G ^ 7 ^ 0 } of the 
following system 


(4.3) 



^ ^ / «! 

dt viieii" 



Pi /3i7 \ 


rj 


is global. Obviously, ri is the explosion time of (4.3), too. Suppose that 
Ti < 00 . On [0,ri), we put X = Then, it is 

easy to verify that {X,Y,Z) satisfies X{t) > 0,Y{t) > X‘^{t)Z‘^{t) and also 
satisfies the following equations 


(4.4) 


' dX 
dt 

dt 

dZ 

> dt 


= -X^Z, 

= -AaiXPZ + Aai-fX^Z - A{PiXP + /3ijX‘i)Y, 

= Y - 2aiXP-^ + 2ai7X9-2 - 2{PiXP + /3i7X'?)Z. 




























6 


TAKESHI UCHITANE, TA VIET TON AND ATSUSHI YAGI 


Furthermore, 

(4.5) limsup[X(t) + y(t) + \Z{t)\ + = oo. 

4—>-ri 

By introducing a function 

H = + X'i-^ + + y + y-4 + m 

with a sufficiently large M > 0, we observe that 

rfzj 

+ 2Y[-AaiXPZ + Aai-fX’iZ - 4(/3iXP + /3i7X'')y] 

+ 2MZ[Y - 2aiXP-'^ + 2ai-iX‘i-'^ - 2(/3iXP + ;8i7X'')Z] 

- 4aiXPZ + 4ai7y''Z - 4(/3iyP + (5i-iX^)Y + 4X-2Z 

= -{q- 2)X'^Z - SaiXPYZ + 8ai7y5yZ + 2MYZ - AMuiXP-'^Z 
+ (4Mai7 -q + 4)y‘?"^Z - AaiXPZ + 4ai7X‘'Z 

- 8{PiXP + /3i7X«)y2 - 4M{PiXP + (3ijX'J)Z^ 

- AiPiXP + /3i7X9)y + 4 x-2z. 

It is easily seen that, for a sufficient small e > 0, it holds true that 

eX'^y + e-^X^-'^ > eX'^+^Z^ + > 2e-^X'?|Z|, 

gXPy2 + MPijXPZ"^ > 2^/eMI3i-fXPY\Zl 
gX-Jy^ + M/3 i 7X‘'Z2 > 2yeM/3i7X‘'y|Z|, 

X'?|Z| + z2 + ^ > (X5 + e-^)|Z| >max{X'?|Z|,XP|Z|,X''-2|Z|}, 
X-4 + Z^ > 2X-2|y|,Z^ + > 2M|Z|, 

(X'? + e“^)|Z| > XP-2|Z| (ifp>2), 

(M^X"^ + 1)|Z| > MXP"2|y| (ifp<2). 

In addition, it is clear that + MZ^ > 2MYZ. Then it follows that 

there exists Mi > 0 such that for X > 0, X > X^Z^, Z G M, ^ is estimated 
by ^ < MiH on [0,ri). Therefore, by the comparison theorem, we obtain 

H{t) = X''-^(t)+X''-2(t)+y2(t)+MZ2(t)+y(t)+X-^(t)+M < 


for all t € [0, ri). Thus, due to (4.5), ri = oo. Therefore, the solution of (4.1) 
must be global. □ 
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4.2. Stochastic case (ai + (T 2 > 0). In this subsection, we consider the 
stochastic case. The system (2.2) becomes 


(4.6) 


' dxi = Vidt + aidwi{t), 

I r 


dvi = < - 


ai7 


Xi — Xj\\P 


\xi — Xj\Yi 


/3i 


\Xi - Xj\\P 


+ 


/3i7 


\Xi — Xj\[i 


{Xi - Xj) 

(vi - Vj)'^dt, 


where i,j = l,2,i ^ j. For (4.6) the situation is not similar to that of the 
deterministic case. Precisely, if d > max{g — 4, 2} and q > 2 then the global 
existence is shown, while if d = 1 or 2 then some solution may explode at a 
finite time. 

Theorem 4.2. Let d > max {(7 — 4,2} and q > 2. Then, for any initial 
value (x^,v^) G M(2), (4.6) has a unique global solution in ]R(2). 

Proof. From Theorem 3.1, there exists a local solution of (4.6) defined 
on [0, r}), where r} is an explosion time. In that interval we have 

' Vi{t) + V2{t) = Ui(0) + U2(0), 

Xi{t) + X2{t) = [ui(0) + U2(0)]t + aiWi{t) + (72W2{t) + a;i(0) + X2(0), 

d{xi - X 2 ) = (ui - V 2 )dt + aidwi{t) - a 2 dw 2 {t), 


d{vi - U2) = I 


- U 2 ) = I - 2 
-2 


ai 


ai7 


||xi - X2\\P 

/3i 


+ 


\xi — X 2 IP 

hi 


(xi - X2) 


{vi - U2)|dt. 


\\xi-X2\\P ||a:i - X2||''_ 

Then r} becomes an explosion time of the following system 


(4.7) 


d( = Tpdt + adw{t), 

ai7 


dfj = 


-2 


C-2 


h , hi 


+ 


IICIP IICII" 


if 


dt, 


too, where C, = xi — X 2 ,'if = vi — V 2 ,cr = + and w{t) = ^[(XiWi{t) — 

(J 2 W 2 {t)\ is also a d- dimensional Brownian motion in (ll,d^, {J}}t>o,lP). By 
putting 


(4.8) 


X = 


1 


y = 


z = {Ch) 
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and using the Ito formula, it is easily obtained that on [0, r^), {X, Y, Z) with 
X{t) > 0,Y{t) >0 satisfies the equations: 

(4.9) 

' dX =[- X^Z - dt - aX^iC, dw), 

' dY =[-4.aiXPZ + AaijX'^Z - 4 {PiXP + l3i-fX’^)Y]dt, 

, dZ =[Y - 2aiXP-^ + 2ai7X'?-2 _ 2{I3iXP + ^i-fX‘^)Z]dt + a{^l;,dw). 
Let us define a sequence of stopping times by putting, for each integer k > ko, 

,k) or Y{t) ^ [0,fc)| , 

where fco > 0 is a sufficiently large number such that (X(0),y(0)) G (^,/co)x 
[0, feo). We here use convention that the infimum of the empty set is oo. Since 
Tk is nondecreasing as k ^ oo, there exists a limit Too = hmfc_>.oo r^. It is clear 
that Too < Ti a.s. We can in fact show that Too = oo a.s. Suppose the con¬ 
trary, then there would exist T > 0 and s G (0,1) such that P{too <T}>s. 
By denoting = {t^ <T}, there exists ki > ko such that 

(4.10) P(flfc) > e for all k > ki. 

Consider the following function in x M,M+) : 

(4.11) V{X,Y,Z) =X^ + X-'^+ Y‘^ + MZ^ + M, 


Tfc = inf t > 0:X(t) ^ (- 


where M > 0 is a sufficiently large number and 0 is a fixed exponent such 
that max{g — 6,0} < 0 < min{d — 2,q — 2}. If (X{t),Y(t), Z{t)) G x M, 
by using the ltd formula, we get 


(4.12) 


dV{X{t), Y{t), Z{t)) =fiX{t),Y{t), Z{t))dt 

+ {9{X{t), y {t), Z{t),Cit),^p{t)),dw{t)). 


Here, 

(4.13) 

f{X,Y,Z) = 

eX^-^[-X^Z - + 4:X-^Z 

+ 2{d + 2)a^X-^ + 2Y[-AaiXPZ + Aai-fX^Z - 4(/3iXP + /3i7X'')y] 

+ 2MZ[Y - 2aiXP-‘^ + 2ai7X'?"2 - 2{I3iXP + l3i'yX^)Z] + a‘^MY 
= - eX^+'^Z - 4aiMXP-'^Z + 2ai7MX‘?-2^ - SaiXPYZ + 8ai7X‘'yy 


+ 2MYZ + Mcr'^Y + 4X-^Z + 2{d + 2)o-^y ^ 


{d-2- 9)ea^ 


X«+2 _ 8{PiXP + Pi-fXi)Y‘^ - 


AMi^iXP + f3i-fX‘i)Z^. 


2 
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And 5 is a suitable function. As for the deterministic case, it holds true that 
/I 






M/?i7 


> 2eX^—|Z|, 


eAPy2 + M/3 i 7XPZ2 > 2v^eM^APy|Z|, 
gjs^gy2 M^xlX^Z'^ > 2v^eM/3i7A''y|Z|, 

y2 + ^2 > 2Y\Z\,Z‘^ + 1 > 2\Z\,Y‘^ + 1 > 2Y,X-^ + M'^ > 2MX-‘^, 

with a sufficiently small e > 0. Whenp > 2, since '^'*' 2 '*'^ > max{q—2, 0+2} > 
p — 2 > 0, we have 

2eX^^^\Z\ + Mi\Z\ > max{M‘^X^+‘^\Z\,M‘^XP-^\Z\,M^X‘i-‘^\Z\} 
with a sufficiently large Mi > 0. Meanwhile, when 1 < p < 2, we have 
XIX-"^ + 2MZ^ + M^> 2MX-'^\Z\ + 2M^\Z\ > 2MXP-^\Z\. 

Thus, whatever p is, there exists M 2 > 0 such that 

fix, Y, Z) < M 2 V{X, Y, Z) for every A > 0, T > Z G M. 

Since for every t > 0 it holds true from (4.8) and (4.12) that 

{X{t ATk),Y{t ATk), Zit ATk)) G R+ xM, y(t Arfc) > A^(t Arfc)Z^(t Ar^), 


we have 

r-sArj, 


ps/\'rfin ps/\Tf^ 

/ dV{X{t),Y{t),Z{t))< M2V(X{t),Yit),Z{t))dt 

Jo Jo 

+ / {giJ<:{t),y{t),z{t)),dwi{t)). 

Jo 

Taking the expectation on both side of this inequality gives 
(iEy(A(sATfc),y(sATfc),Z(sArfc)) < M 2 Ey(A(sATfc), y(sATfc), Z(sATfc))(is, 
from which it follows that for every s > 0, 

Ey(A(s A Tk),Y{s A Tk), Z{s A r^)) < y(A(0), y(0), Z{Q))e^^+ 

In particular, 

(4.14) Ey(A(T A Tk),Y{T A Tk), Z{T A Tk)) < y(A(0), T(0), Z(0))e'^2'^. 
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On the other hand, for every uj G X{Tk){uj) G {k, or Y{Tk){oj) = k. 
Then, 

V{X{T ATk),Y{T ATk),Z{T Ark)) > ak, 
where ak = min | A:^, . Combining this with (4.10), we obtain that 

EV{X{T ATk),Y{T ATk),Z(T Ark)) > E[ln^V{xTATk,yTATk)] > e^fc- 

Therefore, due to (4.14), V{X{0),Y (0), Z{0))e^^'^ > eak- Letting A: — > oo we 
arrive at a contradiction oo > V{X{0),Y (0), Z{0))e^^'^ > oo. Thus Too = oo 
a.s. and consequently, = oo a.s. The proof is now complete. □ 

5. Numerical examples. In this section, we present some numerical 
results. First, we give examples which shows robustness of hsh schooling; 
second, examples which suggest possibility of collision. 

5.1. Robustness. Let us hrst observe examples that show that, if cjj are 
all sufficiently small, then the schooling is strongly robust. 

Set a = 1, (3 = 0.5, r = l,p = 3,q = 4:,ai = 0.015, and Fi{t, Xi,Vi) = —5vi. 
We consider 100(= N) particles in the d-dimensional space, where d = 2, 3. 
An initial value ®(0) is generated randomly in [0, lO]^*^^*^ and i;(0) = 0. 
Figure 1 illustrates positions of particles and their velocity vectors at t = 
0,5,10,15 in M^. Figure 2 does the same at t = 0,10,20,30 in M^. 
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Initial state t =5 




t =10 


t =15 



Fig 1. Schooling in 2-dim. space 
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Initial state 


t=10 






Fig 2. Schooling in 3-dim. space 


5.2. Collision. Let us next observe examples suggesting collision of two 
particles in the d-dimensional space, where d = 1,2, with sufficiently small 
initial distance when Uj are not so small. 

For the case d = 1, we set a = 5, /3 = 1, r = 0.5, p = 3, q = 4, Ui = u, and 
Fi{t, Xi, Vi) = —Vi. An initial value a;(0) is generated randomly in [0,1]^ and 
ri(0) = 0. Figure 3 illustrates trajectories of two particles when u = 0,0.15, 5. 
If a is small (i.e., a = 0,0.15), collision does not take place. Meanwhile if a 
is large (<7 = 5), we observe that collision takes place. 

For the case d = 2, we set a = 7, /3 = 19, r = l,p = 3, g = 4, cTj = 9, and 
Fj(t, Xi,Vi) = —5uj. An initial value a;(0) is generated randomly in [0, 5]'^ and 
u(0) = 0. Figure 4 illustrates behavior of the distance of the two particles 
xi and X 2 - 
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Fig 3. Collision in 1-dim. space 
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Fig 4 . Collision in 2-dim. space 
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